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Abstract 

This work is dedicated to the simulation of full cycles of the electrical activity of the 
heart and the corresponding body surface potential. The model is based on a realistic torso 
and heart anatomy, including ventricles and atria. One of the specificities of our approach 
is to model the atria as a surface, which is the kind of data typically provided by medical 
imaging for thin volumes. The bidomain equations are considered in their usual formulation 
in the ventricles, and in a surface formulation on the atria. Two ionic models are used: 
the Courtemanche-Ramirez-Nattel model on the atria, and the “Minimal model for human 
Ventricular action potentials” (MV) by Bueno-Orovio, Cherry and Fenton in the ventricles. 
The heart is weakly coupled to the torso by a Robin boundary condition based on a resistor- 
capacitor transmission condition. Various ECGs are simulated in healthy and pathological 
conditions (left and right bundle branch blocks, Bachmann’s bundle block, Wolff-Parkinson- 
White syndrome). To assess the numerical ECGs, we use several qualitative and quantitative 
criteria found in the medical literature. Our simulator can also be used to generate the signals 
measured by a vest of electrodes. This capability is illustrated at the end of the article. 
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Introduction 

An electrocardiogram (ECG) is a recording of the electrical activity of the heart pOl IHH] , 
The standard 12-lead ECG is obtained from 9 electrodes located on the body surface. This 
non-invasive and inexpensive procedure is probably the most used clinical tool for the detec¬ 
tion of cardiac pathologies. The motivation for the modeling and the simulation of ECGs is 
twofold. Eirst, the ECG is a simple output of the complex simulation of the cardiac electrical 
activity: while the latter is difficult to validate without invasive measurements, the former is 
easy to assess by a medical expert. Second, the simulation of ECGs, also known as the forward 
problem of electrocardiography [3H] , can be viewed as a step toward the inverse problem of elec¬ 
trocardiography ISD EH SHIES]. Indeed, the inverse problem can be reformulated as a problem 
of identifying the parameters of the model used to simulate ECGs. Erom this perspective, the 
design of a model based on biophysical principles and able to produce ECGs in healthy and 
pathological condition is an important endeavor. This is the main purpose of this study. 

Many attempts at simulating ECGs can be found in the literature m EH EEl [33]. The 
simulation of 12-lead ECG based on partial differential equations (PDE) - as opposed to cellular 
automata [6l] - appeared during the last decade 13 HU [3ni EH EOj. More recently, a focus on 
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the T-wave was proposed [23 El- Simulations of ECGs of patients with Left Bundle Branch 
Blocks were obtained [52]. The precordials leads of the EGG of a rabbit was approximated 
with a variant of the concept of cardiac vector [35]. Because of the difficulty in imaging and 
modeling the atria, all these studies only consider the ventricles. As a consequence, they cannot 
produce the P-wave of the EGG. On the contrary, papers considering atria usually do not include 
ventricles [37] or do not provide ECGs [221 [36] . 

Electrically, the atria and the ventricles are coupled only in a very small area, the atri¬ 
oventricular node. This explains why these studies can concentrate on the upper and lower 
chambers of the heart only. Nevertheless, to produce full ECGs, or to simulate pathologies like 
Wolff-Parkinson-White syndrome, it is necessary to address both atria and ventricles. Recently, 
a full cycle electrocardiogram was proposed for an idealized geometry [56]. To our knowledge, 
the present paper is the first one to propose a full cycle EGG - including the P, QRS and T 
waves - based on a real anatomy of the heart and on the simulation of a PDE-based biophysical 
model. 

Another limitation of the existing works on EGG simulation is their lack of precise evaluation 
of the results. This question is indeed delicate since a “healthy EGG” is not a unique object, 
and there is no obvious metric to measure the discrepancy between this somehow fuzzy notion 
and the result of a simulation. In this paper, we gather many criteria, found in the medical 
literature, that can be used to assess both qualitatively and quantitatively the numerical ECGs. 
We show that our healthy EGG fulfills almost all the qualitative and quantitative properties of 
real ECGs, which is a significant progress with respect to the state of art. Pathological cases 
are also investigated in order to show the capability of our model to predict the features used 
by medical doctors to detect a disease. 

Here is a brief description of our approach. A standard 3D bidomain model is used for 
the ventricles [HIES] and a recent asymptotic surface-based bidomain model is used for the 
atria mm- Two different ionic models are considered: the Courtemanche-Ramirez-Nattel 
model m on the atria, and the “Minimal model for human Ventricular action potentials” (MV) 
in the ventricles [9]. The coupling between the heart and the body is based on a resistor-capacitor 
coupling condition [5]. 

The outline of this article is as follows. In Section 1, the geometry of the heart used for 
the simulation is described. It is based on a surface region for the atria and a volume re¬ 
gion for the ventricles. Its main characteristics are compared with those of a normal human 
heart. Section 2 deals with the biophysical modeling of the atria and the ventricles and the 
coupling condition with the rest of the body. Section 3 concerns the simulations of the stan¬ 
dard 12-lead EGG. A healthy case is given and validated against numerous criteria used to 
assess real electrocardiograms. Some pathological cases are also studied: left and right bundle 
branch blocks, Bachmann’s bundle block, and the Wolff-Parkinson-White syndrome which is 
a pathology caused by the presence of an abnormal accessory electrical conduction pathway 
between the atria and the ventricles. In the last part of this section, we investigate the impact 
of the ionic models on the ECGs by using in the ventricle and in the atria the phenomenological 
Mitchell-Schaeffer model [15]. In Section 4, we show that our simulator can also produce signals 
that are richer than the standard EGG. As an illustration, we analyse the potential measured 
by a “virtual electrode vest” made of 1216 electrodes. In particular, the correlation between 
the signals of different electrodes is studied. We also give an analysis of the dependence of the 
electrode measures with respect to their positions on the body. 
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Left Ventricle 

Volume (ml) 

Mass (g) 

Mitral (cm) 

Aortic (cm) 

Measures 

53.7 

111.3 

2.9 

2.3 

Referenee 

46± 11 

112 ±27 

2.5 ±0.4 

2.3 ±0.2 

Left Atrium 

Major axis (cm) 

Area (cm^) 

Volume (ml) 

Pulmonary Veins (cm) 

Measures 

4.65 

15.9 

47 

1.14-1.45 

Referenee 

3.4 ±0.6 

17.5 ±2.5 

58 ±34 

1.3 ±0.2 


Table 1: Comparisons of the model dimensions with typical end-systolic values found in the 
literature [la Eg sa [13 Eg [la ED] 

1 Whole heart mesh 

To obtain full cycle ECGs, the first step is to build a whole heart realistic mesh. The ventri¬ 
cles can be easily obtained from medical imaging and meshed in 3D. On the contrary, the atria 
have a very thin wall which makes them difficult to image in 3D. In addition, generating a 3D 
mesh on these very thin volumes would dramatically, and uselessly, increase the computational 
cost. For these reasons, we choose to model the geometry of the atria as a surface. We therefore 
obtain an hybrid mesh, made of tetrahedra in the ventricles and of triangles in the atria. 

The heart model was obtained from an anatomical data set called Zygot^ The 3-matic 
software was used to obtain a surface mesh satisfying the standard quality criteria of a finite 
element mesh, and Yams 1201 to refine the surface mesh. Then, the volume of the two ventricles 
was meshed using Gmsh m- We can see in Figure different views of the whole mesh, which 
contains about 230,600 tetrahedra, 73,500 triangles and 67,300 vertices (note that we checked 
that the EGG simulation was not affected by refining the mesh in the ventricles, increasing the 
number of tetrahedra up to 2,511,400). A simplified mesh of the body (Fig. [^, including the 
lungs and the ribs, was also built from the Zygote data set and the aforementioned software. 
The body mesh contains 408,171 tetrahedra, 89,222 triangles and 85,196 vertices. 

The mechanical deformation of the heart is not taken into account in this work. The 
dimensions of the fixed domain correspond to the end of the systole (small ventricles, large 
atria). Tableshows a comparison of a few dimensions of the geometrical model with standard 
end-systolic values. The following quantities are compared: left ventricle volume and mass, 
mitral and aortic valves diameters, left atrium major axis, area, volume and four pulmonary 
veins diameters. We observe a good agreement with the values found in the literature Ha [23112]. 
We also have a good agreement for the diameters of mitral m and aortic [Ml valves, the surface 
of the left atrium [HEg. 

Cardiac tissue has a fiber architecture. The electrical conductivity is higher along the fibers 
than in the transverse direction. This implies that the fiber orientation is very important in 
the study of the electrical activity of the heart. To identify and to prescribe the fibers at 
the endocardium and at the epicardium of the atria, we used [2311211 ESj. As we can see in 
Figure]^ (top), the fibers orientation may vary extremely quickly across the thickness. The colors 
represent the angle 9 defined as half of the angular difference between the endocardium and the 
epicardium. We used [331311 to prescribe the fibers in the ventricles, see Figure(bottom). 

Figure [g represents a schematic view of the heart conduction system in a healthy heart: the 
sinus and atrioventricular nodes, the Bachmann bundle and the Purkinje fibers. In this work, 
the atrio-ventricular node and the Purkinje fibers are not explicitly modeled (see below). 

^www.Sdscience.com 
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Figure 1: Whole heart mesh (right) and body mesh (left). 

2 Modeling assumptions 

In this section, we present the electrophysiology equations and the ionic models used in the 
ventricles and the atria. We also present the coupling conditions between the atria and the 
ventricles and between the heart and the body. 

2.1 Bidomain model 

In order to describe the electrical potential in the heart, we use the standard nonlinear 
reaction-diffusion equations known as the bidomain equations |55^ l58] . In terms of extracellular 
potential Ue and transmembrane potential Vm = Ui — Ue^ with Ui the intracellular potential, the 
bidomain model reads 

+ lioniVm, Wl,---, Wnfj - div(<Tj • VVn) = div(<Tj • V«e) + Amlapp, 

div((CTi + ae) • Vwe) = -div(aj • vv«), 

in B X (0,T), where B denotes the 3D domain of interest, is a positive constant denoting 
the ratio of membrane area per unit volume. Cm the membrane capacitance per unit surface, 
lion the ionic current which depends on n ionic variables wi^ ..., Wn and lapp a given applied 
stimulus current. 

We assume that the heart is isolated, so we make the standard assumption that the extra¬ 
cellular current does not ffow through the epicardium: 

(^e • V?ie) • = 0, ind)Bx(0,T). (2) 


dVm 
dt 


A In ^ 
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Figure 2: Fibers directions at the atrial endocardium (top-left) and atrial epicardium (top- 
right), and in ventricles (bottom). 



Figure 3: Heart conduction system 
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The second boundary condition conies from the fact that, by definition, the intra-cellular current 
does not propagate outside the heart iU: 

(di-Vue)'n = in5;Bx(0,r). (3) 

In order to define the lion term, equations Q must be coupled with a ionic model, i.e. a 
system of nonlinear ordinary differential equations (ODEs). For the ventricular domain, we use 
the MV model, which is a phenomenological model associated with three ionic currents, three 
gate variables, and governed by 28 parameters [9]. 

In order to include the anisotropy between the orthogonal and the tangent directions of the 
fibers, the conductivity tensors Si and Se are defined by 


where / denotes the 3D identity matrix, the vector r is of unit length and parallel to the local 
fiber direction, and and are respectively the conductivity coefficients in the intra- and 
extra-cellular ventricular medium measured along and across the fiber direction. 

The bidomain model can be rewritten in weak form as follows. For all t > 0, find G 

Ue{'^ t) G H^{B) and t),..., Wn{'^ t) G L^{B) with — 0, such that 


(^Cni T lioniy^m-) 5 • • • 5 0 T 


/ {Si + Se) 


V Up 


•V^+ / 
Jb 


Si • VI4 


— I 

JB 

• = 0 , 


+ V?ie) 


• V(/) 


(4) 


for all 0, '0 G H^{B) such that '0 = 0. Under some regularity assumptions, we have existence 
and uniqueness of a solution of the bidomain model The hypothesis fj^Up = 0 is necessary 

in order to have uniqueness and we show in Section 2^ how to adapt this condition when atria 
and ventricles are coupled. Note that we do not model the blood within the atria and ventricles. 


2.2 Surface bidomain model 

As explained in Section it is more convenient to work with a surface mesh for the atria. 
The electrophysiology model used in the present work for the atria was derived from the vol¬ 
ume bidomain model and set on the midsurface. It was obtained from a rigorous asymptotic 
analysis im and was specifically designed for thin cardiac structures [15]. Compared to its 3D 
counterpart, it is very attractive from a computational standpoint. It is worth emphasizing 
that, even if it is set on the midsurface, it actually takes into account volume effects, like the 
strong anisotropy variations across the thickness. In m, we have also performed a numerical 
assessment of the surface model by comparing the results given by the 3D bidomain model and 
the surface model for different geometries and we have used realistic values characteristic of 
atrial electrophysiology for all dimensions and parameters. In particular, we have showed that 
with a thickness of 0.2mm - which can be seen as a standard atrial wall thickness, see [3] - the 
agreement between 3D and surface models is excellent. 

Let S denote the midsurface of the wall and H^{S) the associated functional space. The 
surface-based bidomain model can be rewritten in weak form as follows: for all t > 0, find 
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G i7^(5), Ue{'^ t) G H^{S) and wi{’^ t),..., Wn{'->t) G L'^{S) with Ue = 0, such that 



for all 0, '0 G H^{S) such that '0 = 0. We denote by and the conductivity coefficients 
in the intra- and extra-cellular atrial medium measured along and across the fiber direction. 
The intra- and extra-cellular diffusion tensors and ag are defined by 


7 " I / Q 

gi,e = 1+ K 


a,I 


a,t 

^ie 


) [Ioi0)ljo ® To + JQiG)Zo ® 


( 6 ) 


where / denotes the identity tensor in the tangential plane, Tq is a unit vector parallel to the local 
fiber direction on the atria midsurface, and such that (tq,t^) gives an orthonormal basis 
of the tangential plane. We use the fibers direction at the endocardium and at the epicardium 
to define the fibers direction Tq on the atria midsurface and the angle variation 9 between 
the endocardium and the epicardium. The effect of angular variations appears in the model 
with the coefficients Io{9) = ^ ^ sin(20) and Jo{9) = 1 — Io{9). Note that Jo{9) = 0 (and 

^o(^) = 1) if si^nd only if 0 = 0, which corresponds to a constant direction in the thickness and 
then ai^e — ~ ® To- By contrast, important angular variations make Iq 

decrease and Jq increase in Q and the diffusion becomes more isotropic. This model has been 
compared m to several 3D models proposed in the literature [IH1E21I13]. 

The Courtemanche-Ramirez-Nattel ionic model m is considered. It includes 12 ionic cur¬ 
rents and 20 other variables. The two atria are connected only by two regions, the Bachmann 
bundle and the Fossa Ovalis. We refer to m for more details. 


2.3 Simulations on the whole heart 


Coupled model As seen in Sections 2.1 and 2.2, the unique solution G H^{B) of Q is 

s.t. f^Ue = 0 and the unique solution G H^{S) of is s.t. Tig = 0. Here we consider 

the whole domain BuS and a new global criterion. The resulting coupled problem is well-posed 
at the discrete level, but its mathematical analysis remains to be done. Let = B^US^, where 
Bh is the mesh of the ventricles and is the mesh of the atria, and let be the line such that 
BhFSh — F}i. We denote by the restriction of a function to a subdomain D. The finite 
dimensional approximation space Vh is then defined by: Uh G Vh if and only if Uh is continuous 
in Qh, iBhUh e ISh'^h e and Uh = 0. Using Q, Q, the full model reads, 

find (ue./i., Vm,h) e Vft, such that e Vh, 
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body 

arj. 

^bones 

U rp 

lungs 

U p 

3.010-4 

1.210-4 

2.010-^ 


Table 2: Torso conductivity parameters (all in S.cm 



Figure 4: Left - Surface connection between the atria and the ventricles, Right - Kent bundle 


Connection surface As previously mentioned, the atrioventricular node is the only pathway 
for the electrical signal between the atria and the ventricles. From a physiological point of view, 
a fibrous skeleton separates atria boundaries from ventricles epicardium. This layer isolates the 
atrial cells from the ventricular ones [ 1211111 . We propose to model this fibrous skeleton with a 
thin layer of the atrial surface, represented on the left of Figure]^ The idea is that in this area 
there is only a low conduction of the extracellular potential. In this region, denoted by Sc C 5, 
the intracellular conductivity is set to zero and the extracellular conductivity denoted by 
is very low (see Table [^. Finally, the surface and volume bidomain equations 0 are solved 
simultaneously on this “hybrid” domain with 


CT', 


V.t J . / V. 

I + 


V^l 


V,t\ 




gi,e = - ^ii) [lo{0)rQ (g) To + M0)t^ ® , in 5 \ 5c, 


ai = 0, and ac = in 5c. 


Parameters and applied currents The values of the membrane parameters are = 
200.0 cm“^ and Cm = 10“^ mF.cm“^ for the whole heart. The conductivity takes different 
values depending on the region in the ventricles and atria (Table |^. 

In the atria, the regions of fast conduction are the Bachmann bundle (BB), see Figure 
the Crista Terminalis (CT) and the pectinate muscles (PM). By contrast, the Fossa Ovalis 
(FO) is a region of slow conduction. In order to model the different propagation velocities, we 
modify the values of the maximal conductance of the Na‘^^ current Table gives the 
parameters used for gjsfa- Furthermore, to reduced the action potential duration, the parameter 
gk^ is chosen five times as big as in the original model HZ!. 







v,t 

V,t 

V.t 

v,L 

a,t 

a,l 

a,t 

ad 

S,t 

(Te 

(Te 



CTe 

CTe 



(Te 

4.810“^ 

1.610“^ 

1.610“^ 

1.610"^ 

9.010-^ 

2.510"^ 

2.510"^ 

2.510"^ 

7.510-^ 


Table 3: Conductivity parameters (all in S.cm 


regular tissue 

PM 

CT 

BB 

FO 

7.8 

11.7 

31.2 

46.8 

3.9 


Table 4: Maximal conductance g^a in the different atrial areas (all in nS.pF 


In the ventricles, we modify the duration of the plateau too. In the MV model, we change 
the values of Tgoi parameter in order to reduce the action potential duration for epicardial, 
endocardial and midmyocardial cells. This heterogeneity is considered in the left ventricle, for 
the positivity of T wave |65^ I34|. In the right ventricle, the cells are considered homogeneous 
and their parameters are taken as in the left ventricle epicardium, except for Tsoi (Table [^. 

Activation is initiated at the sinus node with a stimulus of 2ms which triggers a depolar¬ 
ization wavefront in the atria (Figure |^. For the sake of simplicity, the atrioventricular node, 
which is the only electrical connection between the atria and the ventricles, is not modeled with 
a sophisticated physiological model. Instead, the excitation is triggered in the ventricle after a 
parameterized delay (in healthy condition, we choose to start it at 190 ms). Similarly, the fast 
conduction in the Purkinje fibers (Figure]^ is modeled with a predefined stimulus pattern: a 
time-dependent thin subendocardial layer is activated by an external current for 5 ms on both 
right and left ventricles [5]. Figure shows the propagation of the activation in the atria and 
Figure]^ shows the early stage of the activation in the ventricles. 

Simulation results The various simulations of this article are performed with the finite el¬ 
ement library FELiSc^^^ developed at Inria. Problem Q is solved with a BDF2 semi-implicit 
scheme [3]. Figure]^ shows a full cardiac cycle. The corresponding first lead electrocardiogram 
is also represented. The electrical signal starts at the sinus node where the atrial depolar¬ 
ization (AD) begins. By 50 ms the wave finishes spreading along the Crista Terminalis as a 
consequence of the high conductivity in this part, see Figure Importantly, because of the 
rapid conduction in the Bachmann bundle, the wave spreads to the left atrial appendage and 
activates a substantial part of the left atrial wall. The depolarization of the right and left atria 
terminates at 100 ms and 110 ms, respectively. The ventricular depolarization begins at 190 ms. 
During this period, the atrial repolarization (AR) occurs. As we can see in Figure]^ at 200ms 
the endocardium of the ventricles rapidly depolarizes. Then, the wave propagates across the 
ventricles. The repolarization ends at 400 ms in the right ventricle and at 430 ms in the left 
ventricle. 

^http://felisce.gforge.inria.fr 



EPI 

ENDO 

M 

RV 

m 

30.0181 

40.0 

91.0 

/ 

heart 

19.0 

30.0 

45.0 

20.0 


Table 5: Changed ionic parameter Tgoi of MV compared to [9] 
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Figure 6: Action potential in the ventricles. A stimulus is applied from t=190 ms to t=195 ms. 

2.4 Coupling with the body 

The last step in order to obtain an electrocardiogram is to couple the heart model with a 
diffusion problem in the rest of the body 

— div((j7"— 0, in (8) 

where the electrical conductivity ar takes different scalar values in the ribs and the lungs (see 
[To] and Table . 

On the body surface an homogeneous Neumann boundary condition is imposed 

aT'VuT • n = 0. To define the transmission conditions at the heart-body interface we 

assume that the extracellular current does not ffow through the pericardium (isolated heart 
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AD - 100ms AR/VD - 200ms 



Figure 7: Simulations of heart depolarization in a healthy case with the corresponding electro¬ 
cardiogram first lead 
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Figure 8: Body potential simulations of heart depolarization in a healthy case (Figure]^, with 
coupling condition (10) 
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assumption): 

Se • Vl^e • 'u = 0 (9) 

and we consider the resistor-capacitor conditions [5] : 

RpiarVur) ■ n = + {u, - ut) (10) 

where Cp and Rp stand for the capacitance and resistance of the pericardium, respectively. 

Condition is an approximation that is known i-^supposed not to affect too much the shape 
of the ECG, as shown in [5]. It allows us to solve the heart-body system as a one-way coupled 
problem, which dramatically reduces its computational cost. 

Condition is a generalization of the condition that is usually adopted (continuity of 
the potential). It allows us to model the fact that the transmission of potential through the 
pericardium is not perfect, and can be different for the ventricles and the atria. We take 
Rp = 10^ fi.cm^ on the surface in contact with the ventricles and Rp = 10^ ri.cm^ on the 
surface in contact with the atria. These values are chosen to obtain correct relative amplitudes 
of the P and R waves. We neglect the capacitor effect by taking Cp = 0 mF.cm^ in (10). 


The transmission between the heart and the body is therefore modeled as a Robin boundary 
condition: 

Rp{o't^Ut) • fi -\- Ut — (11) 

Conditions (§ and ( [Tol ) are further commented in Section]^ Figureshows the body surface 
potential corresponding to the simulation shown in Figure!^ 


2.5 Electrocardiogram computation 


A standard electrocardiogram is based on the body surface potential recorded by 9 elec¬ 
trodes (r ecg — {R^ F, Vi,..., Ve}? see Figure [^. These measures are combined to define 12 
differences of potential, known as the 12 leads of the standard ECG: 


I = ut{L^ ~ ut{R) 

II = ut{F) — ut{R) 

III = ut{F) - ut{L) 
VI = Ut{Vi) - Uyj 

V2 — ut{V2) — 

V3 = uriVs) - 


aVR — 1.3 {ut{R) — Uyo) 
aVL = l.b{uT{L) — Uyo) 
aVF = 1.5{ut{F) — u^) 
V4 = ut(V 4) — Uyj 
F5 = ut{V^) — u^ 

V6 = ut{Vq) - Uy, 


where Uyj — ^{ut{L) -h ut{R) + ut{F)) is the Wilson potential [iOj . 


3 Numerical healthy and pathological ECGs 

In this Section, we present the ECGs provided by the aforementioned model in healthy and 
pathological conditions. The healthy ECG is obtained by carefully choosing the parameters of 
the model in order to match most of the qualitative and quantitative features of a physiological 
ECG. To obtain the pathological ECGs, the approach is different: starting from the nominal 
values corresponding to a healthy ECG, we modify the parameters in order to model the physical 
characteristics of the pathology. Then we observe the effects of these modifications on the 
numerical ECG, and we compare its features with the ones described in the literature. It is 
important to emphasize that, for the pathological cases, the parameters are not intentionally 
fixed to match a given ECG. Thus, if the ECGs obtained after modeling the diseases match the 
main features observed on real patients, it gives confidence in the prediction capabilities of the 
model. 
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I = ut{L) — ut{R) 



Figure 9: Standard 9 electrodes locations and first and second ECG leads. 


3.1 Healthy electrocardiograms 

Fignre [TT] shows the simulated electrocardiogram in healthy conditions, corresponding to the 
simulation of Figure An electrocardiogram is typically described by distinguishing five events 
during the heartbeat, called P, Q, R, S and T “waves” (we will keep this standard terminology 
even though these events have nothing to do with waves). The P wave corresponds to the 
atrial depolarization, the QRS complex corresponds to the ventricular depolarization, the T 
wave corresponds to the ventricular repolarization. The typical durations of each wave, or each 
interval, are given in Table 

Table also presents the durations of the simulated healthy ECG of Figure These 
durations are obtained in the numerical ECG from the landmarks defined according to the 
following rules: the P wave (resp. QRS complex) starts if 1% of the atria (resp. ventricles) is 
activated (i.e. if the transmembrane potential Vm is greater than a threshold voltage Vgate)] fhe 
P wave (resp. QRS complex) ends when 99% of the atria (resp. ventricles) are activated; the T 
wave starts when 20% of the ventricles are repolarized {i.e. Vm < 0); the T wave ends when 99% 
of the ventricles are fully repolarized {Vm < 0). If the minimal value of Vm is Vmin — —80mV and 
its maximal value Vmax = 20mV, we define Vgate = —67mV, which corresponds to a threshold 
voltage 6yo — 0.13 in the MV model. We remark that the waves durations are in general in good 
accordance with the values found in the literature, even though the QRS complex is slightly 
shorter than expected. 

Table gives the main features of each wave in a normal electrocardiogram. Note that the 
simulated ECG fulfills almost all the expected criteria. We only observe a discrepancy in the 
aVL lead, but this lead is not the most important one for the ECG interpretation. 

To qualitatively assess the waves amplitude and orientation, the schematic presented in 
Figure 12 is very convenient. It is adapted from [63] and shows the normal variations of wave 


amplitude measured in adults. A visual comparison of Figures [TT] and [T^ shows that, for almost 
every lead, each wave of our numerical ECG is in the range of the normal values. Note that in 
Figure the length of each wave was arbitrarily chosen as its maximal normal duration. This 
is the reason why the full PQRST duration is so long in this schematic. 

Here is another qualitative assessment. The R wave is known to have an important property 
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P 

wave 

PR 

interval 

Q 

wave 

QR 

interval 

s 

wave 

QRS 

interval 

QT 

interval 

Typical 

Duration 

< 0.12 

0.12 

to 0.21 

< 0.04 

< 0.03 VI - V2 

< 0.05 y5 - V6 

< 0.04 

< 0.10 

0.35 

to 0.45 

Healthy 

ECG 

0.08 

0.19 

0.015 

0.025 VI - V2 

0.02 P5 - V6 

0.01 

0.04 

0.25 


Table 6: Durations of the simulated healthy ECG of Figure [T^ compared with typical durations 
[63] (all in s) 


Wave/Interval 

Description 


Simulated ECG 


< 0.25mV 


/ 0.2mV 

P wave 

positive I, II, V3 to V6 


/ 


negative aVR 


/ 


limb leads < 25% of R 


/ 

Q wave 

precordial leads < 15% of R 


/ 


always negative 


/ except for aVL 


limb leads < 2mV 


/ 

R wave 

precordial leads < 3mV 


/ 


always positive, negative in aVR 

/ 


R wave progression, see Figure 

10 

/ 


always negative 


/ 

S wave 

small I, II, V5, V6 


/ 


important VI to V3 


/ 


—0.05mV to O.lmV 


/ 

ST interval 

isoelectric 


/ 


displacement of 0.02mV in VI, 

V3 

/ 

T wave 

positive I, II, V3 to V6 


/ 


negative aVR (follow the QRS) 

/ 


Table 7: Criteria for a typical electrocardiogram [63] compared with simulated ECG of Figure [TT| 


in the precordial leads: it uniformly progresses from a rS complex (small R wave, deep S 
wave) in VI-V2 to a QRS complex in V5-V6 via a RS complex in V3-V4. This represents the 
rotation of the cardiac vector from right to left, and from back to front and again back. The 
idealized behavior is represented in the top of Figure [To| extracted from [63|. The bottom of the 
same Figure shows the results of our simulation. Given the difficulty to reproduce this subtle 
dynamics with the simulation of a complex biophysical model, the agreement can be considered 
as satisfactory. 

A last qualitative comment is in order. We note that the P wave presents some oscillations 
in all the leads of the numerical ECG. The explanation of these oscillations is the brutal changes 
of the fibers’ direction in the atria. It is also possible that the surface representation of the atria 
accentuates these oscillations by diminishing the natural diffusion arising during the propagation 
of the depolarization front through the thickness. 

3.2 Pathological electrocardiograms 

In this Section, we modify the protocol of the simulation that provided the healthy ECG 
(Figure 0 in order to simulate different cardiac pathologies. Then we verify if the numerical 
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Figure 10: R wave progression in the precordial leads: schematic view from [63] in the top, and 
simulated ECG in the bottom 



Figure 11: Healthy electrocardiogram corresponding to simulation of Figure (voltages (mV) 
versus time (ms)) 
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Figure 12: Representation of the variability of the amplitude observed in adults healthy ECGs, 
according to [63] . 


ECGs present the main features that allow a medical doctor to detect the pathology. The 
different pathologies are schematically represented in Figure [T^ along with the most important 
leads in each case. 


3.2.1 Left and right ventricular block 


We start with a left or a right bundle branch block. In a healthy case, the right and the 
left ventricles are activated simultaneously. Now, in order to simulate a left (or a right) bundle 
branch block, the initial activation is blocked in the left (resp. right) ventricle. On the left of 


Figure 13, we can see a left ventricular block. In order to obtain a left (resp. right) bundle 
branch block, the depolarization of the left (resp. right) Purkinje fibers is delayed [5j. Results 
are reported in Figure for the left and right bundle branch blocks. We recognize the main 
characteristics reported in the medical literature: larger QRS, lead VI without Q-wave m, 
leads VI and V6 similar to those presented in Figure [T^ (left). The QRS-complex exceeds 0.1 
seconds in both cases. Furthermore, it can be seen in Figure [14| that the duration between the 
beginning of the QRS complex and its last positive wave in VI (resp. V6) exceeds 0.04 seconds 
which is a known sign of right (resp. left) bundle branch block [40] . 


3.2.2 Bachmann’s bundle block 

In the heart conduction system, the Bachmann bundle connects the left atrium with the 
right atrium and is the preferential path for the electrical activation of the left atrium. A block 
of the Bachmann bundle is represented in the middle of Figure It is characterized by the 
presence of P-wave duration that equals or exceeds 0.12 seconds and presents usually a bimodal 
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Left ventricular block Bachmann’s bundle block Wolff-Parkinson-White syndrome 



Figure 13: Different pathologies 


morphology, especially in leads I, II, aVF and the lead III becomes biphasic, as we can see 
in Figure 13, This is a very specific sign of left atrial enlargement @01117]. We simulate it 
by decreasing the maximal conductance = 7.8 in the Bachmann bundle. The results are 
given in Figure!^ The more important the block, the more negative the P wave on lead III. 


A negative P wave in the third lead corresponds to the retrograde depolarization of the left 
atrium. The morphology of the simulated P wave is in very good agreement with the criteria 
given in the literature for various degrees of Bachmann’s bundle blocks [HE]. 


3.2.3 WolfF-Parkinson-White syndrome 

The Wolff-Parkinson-White syndrome is one of the numerous pathologies of the conduction 
system of the heart. It corresponds to a pre-excitation syndrome and is caused by the presence of 
one or more abnormal electrical conduction pathways between the atria and the ventricles, called 
bundles of Kent. Electrical signals travel down these abnormal pathways and may stimulate the 
ventricles prematurely. In the right of Figure [T^ we can see a schematic of the Wolff-Parkinson- 
White syndrome. Here we consider one Kent bundle only and we model this abnormal pathway 
by stimulating a ventricle area near of the atria represented at the right of Figure The 
Wolff-Parkinson-White syndrome is commonly diagnosed with the electrocardiogram [Ml. It 
is characterized by a delta wave, a slurring of the initial segment of the QRS complex, due to 
the arrival of the impulses in the ventricles via the abnormal route, which is associated with a 
short PR interval. Another feature is a QRS complex widening with a total duration greater 
than 0.12 seconds. We can indeed observe these characteristics, in particular the delta wave in 
Figure [T6| 

3.3 Comparison with the Mitchell-Schaeffer model 

In this section we are interested in the impact of the ionic model on the ECG simulation. 
The membrane current is now described with the Mitchell-Schaeffer model |45j which is a one- 
current phenomenological ionic model, offering interesting properties with a very limited number 
of parameters. 
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(b) Right ventricular bundle brunch block. 


Figure 14: Left and Right Bundle Brunch Block - Healthy case in red, LBBB in blue and RBBB 
in green (voltages (mV) versus time (ms)) 
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Figure 15: Bachmann’s Bundle Block - Healthy case in red, BBB in blue (voltages (mV) versus 
time (ms)) 



II III 




Figure 16: Wolff-Parkinson-White syndrome - Healthy case in red and WPW in blue (voltages 
(mV) versus time (ms)) 


The Mitchell-Schaeffer model is applied with the same conductivity parameters (except for 
some atrial areas, see below) and the same initial stimulus as described above. Table gives 
the value of the Mitchell-Schaeffer parameters, appropriately rescaled [5]. In order to correctly 
reproduce the T wave, we take into account three layers of cells in the left ventricles and an 
homogeneous tissue in the right ventricle, as described in Section |2.3[ The Tdose parameter 
varies according to the type of cell (Table [^. In the atria, the repolarization “propagates” in 
the same direction as the depolarization. We therefore take a constant value for Tdose^ equal to 
100ms. As previously explained, we changed the values of the maximal conductance g^a in fhe 
different atrial areas in the Courtemanche-Ramirez-Nattel model in order to take into account 
these bundles. The Mitchell-Schaeffer model does not allow the same ffexibility, so we decide 
to directly modify the value of conductivity parameters. The atrial conductivities are modified 
as reported in Table in order to represent the different slow and fast bundles. 


'T~in 

(cm^ 

'^out 

.mA“^) 

'^open 

^endo 
' close 

Mcell 
' close 

(ms) 

epi 

' close 

^RV 
' close 

^min ^max 

(mV) 

^gate 

4.0 

90.0 

300.0 

120.0 

100.0 

80.0 

90.0 

-80.0 20. 

-67.0 


Table 8: Mitchell and Schaeffer parameters and constants (different values of Tdose nre given 
because of an heterogeneous tissue is considered 151) 

Figure [T7| shows the ECGs obtained with the Mitchell-Schaeffer (MS) model and the com¬ 
bined Courtemanche/MV model. We can see that the results are reasonably close. With the 
MS model some oscillations in the P wave and the QRS complex are fixed, but the R wave 
progression in precordial leads is less satisfactory and the T wave in V2 and V3 are not correct. 
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regular tissue 

PM 

CT 

BB 

FO 

a2 

a.t 

(Te 

ad 

ad 

(de 

2.510-^ 

9.010“^ 

2.510“^ 

2.510“^ 

4.510"^ 

1.3510-3 

4.510-3 

4.510-3 

7.510-^ 

2.710-3 

1.0910-2 

1.0910-2 

1.1910-3 
4.310-3 
1.8610-2 
1.8610-2 

2.510-^ 
9.010-^ 
2.2710-3 
2.2710-3 


Table 9: Atrial conductivity parameters (all in S.cm for the Mitchell-Schaeffer model 


In spite of these flaws, those results are reasonable, and could probably be improved with a finer 
choice of the parameters. It is interesting to note that the results of the simulations are robust 
with respect to the choice of the ionic model: the Courtemanche/MV model in general gives 
better results, but it can be replaced by the MS model in order to reduce the computational 
costs and the model complexity without affecting too much the ECG. This remark is especially 
important if the ECG simulator has to be used for inverse problems: in that case, a model with 
a reduced number of parameters may be more attractive. 



II 



III 



aVR 



aVL 



aVF 




Figure 17: ECGs obtained with different ionic models - Courtemanche/MV ionic model in red 
and Mitchell and Schaeffer model in blue (voltages (mV) versus time (ms)) 


4 Virtual “electrode vest” 

In the previous sections, we focused on the 12-lead ECG because it is widely used in practice, 
and it can be easily assessed with the medical literature. But our simulator can of course provide 
more sophisticated measurements, like those obtained with electrode vests. 
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Electrodes vest 2 clusters 3 clusters 



Figure 18: Human body mesh, “electrode vest” (left) and clustering results (right). From left-top 
to right-bottom clustering agglomeration, from 2 to 6 clusters. 


Many studies have been carried out on this topic: on the forward problem and the analysis 
of the number of electrodes [28^ 1^ . or on the reconstruction of the potential on the heart 
surface [S31EH]. Our objective here is less ambitious: we just show an example of a statistical 
analysis that can be done with the measures provided by our EGG simulator. 

To do so, we simulate a virtual “electrode vest” which contains Necg — 1,216 electrodes. 
Figure [T^ (left) shows the “electrode” locations, which correspond to all the nodes of the mesh 
in the red region of the torso. The heart geometry used in this study contains = 28, 510 

boundary vertices. We compute the body surface potential as described in Section 

We are interested in analyzing the electrode signals with respect to their positions on the 
body. In order to divide the signals into different groups, an agglomerative hierarchical clus¬ 
tering analysis is applied to the 1,216 measures. Clustering is a statistical technique used to 
classify data based on similarities^ or distances^ between them. In particular, in agglomerative 
hierarchical clustering method, at the beginning each individual (data) represents a group itself. 
Then, these groups are merged together according to their decreasing distance. The procedure 
is schematized in Figure the bottom of the clustering tree represents the 1, 216 signals, the 
top represents a unique group. In order to measure the similarity between data, an Euclidean 
distance is used 

distCVi, Vj) = /E (12) 

y n 

where 1^, i = 1, • • •, 1216, represents an electrode signal. The criterion used to merge groups 
{linkage) is based on the distance between two groups, defined as the maximal distance between 
the data of the first group and the data of the second one 

dist(C'i, Cj) = max (dist({14 G O}, {I 4 G Cj})) (13) 

k,h 

where C* represents the i—th. cluster [32]. 
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Figure 19: Principal components analysis on electrode vest signals 
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Figure 20: Clustering agglomeration of electrodes measures potential using euclidean distance 
and complete linkage. Colors represent the clusters of points as shown in Figure fig:body-mesh- 
cluster. 


In Figure p!8| we show the agglomeration procedure from 2 to 6 clusters. First, the division 
into two clusters indicates a separation between the heart region (purple area of Figure 18) 
and the rest of the vest (green region). Second, a division into 2 parts underlies the atria- 
ventricular axes (separation of green and yellow zones). Then, we can see from the clustering 
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agglomeration plot of Figure 20, that the separations into 4, 5 and 6 groups (blue, orange and 
red areas) are made at very closed distances. The last subdivision we consider is made of 6 
groups. In Figurewe plot the signals belonging to the 6 clusters and their center values. The 
center values of each cluster are computed as the point minimizing the distance between itself 
and the other points of the cluster 


rrii = argmin ( ^ distCt^-, 14)), Vi. 
^ k 


(14) 


Cluster #1 


Cluster #2 Cluster #3 




Figure 21: Signals in different clusters as shown in Figure [T^right). Black lines indicate the 
“mean” point (points of Figure [T^right)) 


Then, we re-apply the principal components decomposition on the 6 centers of the cluster, 

we 


the points indicated in Figure [Tq (right). Comparing results of Figure 22 with Figure 19 


observe that the first principal component is again much larger than the other ones, and the 
same curves represent the first, second and third principal components. 

The two approaches - the principal components decomposition and the hierarchical cluster¬ 
ing analysis - suggest that it is not necessary to have a very high number of skin electrodes to 
describe the body surface potential. A limited number of correctly positioned electrodes should 
be enough to represent most of the features of the signal. A deeper analysis of the number and 


24 

























the locations of the electrode goes beyond the scope of this study. The purpose here was just to 
illustrate that our simulator could be used to generate more general signals than the standard 
12-lead ECGs. 


Principal components Cumulative eigenvalues 



Figure 22: Principal components analysis on 6 electrodes signals. 


5 Limitations 


This work can be considered as a step forward towards the simulation of ECGs with a 
biophysical model. However, it requires further improvement and several limitations should be 
mentioned. 

A strong coupling with the torso should be done to assess the impact of the isolated heart 
assumption It was indeed shown in [5] that this assumption modified the amplitude of the 
QRS with a factor two without affecting the shape significantly. But it is possible that the 
isolated heart assumption has a stronger effect with other activation sequences. 

The role of parameter Rp in (10) is not completely understood. To get correct relative 
amplitudes of the P and R waves, we took ad hoc values, higher in the atria than in the 
ventricle. We postulated that this might be due to a difference of conductivity at the interfaces 
ventricle-torso and atria-torso but we did not find any evidence in the literature to support 
this assumption. This choice should therefore be seen as a research hypothesis which has to be 
further investigated. 

Another limitation that should be mentioned is the rough modeling of the atrio-ventricular 
node and the Purkinje fibers. 

At last, by adjusting some parameters of the models, it is likely that the realism of the 
simulated ECGs could still be improved. This task, which was mainly done manually in this 
work, should be addressed with an inverse problem strategy mm- 
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6 Conclusions 


We have presented a comprehensive model for the simulation of full cycle ECGs. The 
main ingredients are: volume bidomain equations for the ventricular part coupled with the MV 
ionic model; surface bidomain equations for the atrial part coupled with the Courtemanche 
ionic model; one-way coupling between the heart and the torso through a resistor transmission 
condition. 

This model has provided a healthy EGG whose qualities and flaws have been assessed with 
several qualitative and quantitative criteria extracted from the medical literature. Four patho¬ 
logical cases have been investigated: left and right bundle branch blocks, Bachmann’s bundle 
block, and Wolff-Parkinson-White syndrome. In the healthy case, we have shown that similar 
results can be obtained with the one-current phenomenological Mitchell-Schaeffer model. The 
simulation of a virtual vest of electrodes has illustrated a potential application of this simulator. 

As a natural continuation of this study, the proposed electrical model of the heart could also 
be used in electromechanical simulations, as was done in nam, including a mechanical model 
of the atria. 
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